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Abstract We study the influence of particle shape anisotro- 
py on the occurrence of avalanches in sheared granular me- 
dia. We use molecular dynamic simulations to calculate the 
relative movement of two tectonic plates. Our model con- 
siders irregular polygonal particles constituting the material 
within the shear zone. We find that the magnitude of the 
avalanches is approximately independent on particle shape 
and in good agreement with the Gutenberg-Richter law, but 
the aftershock sequences are strongly influenced by the par- 
ticle anisotropy yielding variations on the exponent charac- 
terizing the empirical Omori's law. Our findings enable one 
to identify the presence of anisotropic particles at the macro- 
mechanical level only by observing the avalanche sequences 
of real faults. In addition, we calculate the probability of oc- 
currence of an avalanche for given values of stiffness or fric- 
tional strength and observe also a significant influence of the 
particle anisotropy. 



1 Introduction 

Natural earthquakes are one of the most catastrophic events 
in nature f\\\ with deep social implications, in terms of hu- 
man casualties and economic loss. Considerable efforts have 
been made to understand the earthquake dynamics and the 
underlying mechanisms prior to their occurrence ||2;[l;3;|5l], 
either through experimental studi es laH jal or through par- 
ticle based numerical models |[ 9 ; [iQ : [ill: [l2n of the relative 
motion of tectonic plates Jl]; ITx l4\ . 

However, in most of the existing numerical models of 
earthquake fault the gouge is represented by discs ifTTI : [l2h 
or spheres flSh . The dynamics of such material within the 
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fault is thought to control the stick-slip instability that char- 
acterizes earthquake process. An understanding of its prop- 
erties is, therefore, vital to understand earthquake dynam- 
ics fl^. For instance, the existence of the gouge within the 
fault has been proposed to explain the low dissipation on 
shear zones, providing new insight into the heat flow para- 
dox In this case, the reduction of the macroscopic fric- 
tion and consequently, the heat generation is attributed to the 
deformational patterns such as rolling of particles ifTTI: [l2ll . 
In laboratory experiments by Maron [g], the influence of par- 
ticle characteristics has also been studied. They found that 
frictional strength and stability of the granular shear zone is 
influenced by particle shape, size distribution and their evo- 
lution through particle crushing. Therefore, to model fault 
gouges one must also include different grain characteristics. 

In this paper, we use a model of polygonal particles fisl : 
[19I1 to address the influence of anisotropy in granular media. 
We study the situation of two tectonic plates with bound- 
aries parallel to the direction along which the tectonic plates 
move - so-called transform boundaries fll: [l4ll . One of the 
most well known examples of such boundaries is the San 
Andres Fault in California where the Pacific plate and the 
North American plate are moving in opposite directions. In 
this particular case, the relative motion of the plates is about 
40 mm/year, and thus the strain accumulation rate is around 
3-10~^ per year ll20ll . Further, different from previous studies 
iflOll . our model considers anisotropic particle shapes. The 
response of the system is characterized by discrete events 
or avalanches whose size covers many orders of magnitude, 
similar to the so-called crackling noise of physical systems 

M. 



We find that the magnitude of the avalanches is indepen- 
dent of the particle ^hape and in agreement with the Guten- 
berg-Richter law [22;]. On the contrary, the distribution of 
the waiting times described by the Omori's law ll23ll strongly 
depends on shape anisotropy. From this result, we raise the 
hypothesis of identifying at the macro-mechanical level the 
presence of anisotropic particles within the gouge, only by 
studying the temporal avalanche sequences. We further ar- 
gue that the existence of this anisotropic gouge in fault zones 
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Fig. 1 Schematic representation of a particle contact. Tiie overlapping 
area a is indicated by the shaded zone. 

might also explain the variation of the decay of the after- 
shock sequences observed in nature. 

In addition, we also compute the conditional probability 
for an avalanche to occur, and found that it decreases log- 
arithmically with the stiffness. This exponential decay also 
depends on particle shape, since anisotropic samples are able 
to mobilize a higher frictional strength when compared to 
isotropic samples. For a given value of mobilized strength 
anisotropic samples also exhibit lower probability of failure. 
Finally, we propose some microstructural features that could 
help to explain the occurrence of avalanches. 

We start in Sec.|2]by describing in detail our model of 
irregular particles as well as the details of our numerical ex- 
periment. In Sec.[3]we characterize and study the system re- 
sponse. In Secs.|4]and|5]we address the influence of particle 
anisotropy on the frequency distribution of avalanches and 
on the width of the time interval where aftershocks occur. 
The weakening and stability of the system is investigated in 
Sec.|6l and in Sec.|7]the main conclusions are discussed. 



2 The model 

We consider a two dimensional model of convex polygons to 
represent the grains of the granular material on a mesoscopic 
scale. The samples consist of isotropic and anisotropic par- 
ticles in order to study the influence of particle shape aniso- 
tropy on the response of a granular packing under very slow 
shear. 

The deformation of the grains is modeled by letting them 
overlap, as sketched in Fig. [T] When two polygons overlap, 
the intersection points between their edges can be defined. 
The segment that connects these points. Pi and P2, gives 
the contact line S — PiP2- The contact force is given by 



f ^ = f + f where f and f" are the elastic and viscous 
contribution. 

The elastic part of the contact force is decomposed as 

= /n"-*^ + ft^'^' where n"^ and are the unitary vectors 
perpendicular and parallel to the contact line S respectively. 
The normal elastic force is calculated as = — fc,i(5, where 
fc„ is the normal stiffness, and 5 the deformation length de- 
fined in terms of the overlapping area a and the length of 
the contact line, 6 ~ a/|S|. The friction force is given by an 
elastic force — —kt^ proportional to the elastic displace- 
ment ^ at each contact, with kt the tangential stiffness. The 
elastic displacement ^ is updated as ^ = ^' + v^At, where ^' 
is the previous length of the spring. At is the time step of the 
molecular dynamic simulation, and Vj the tangential compo- 
nent of the relative velocity at the contact. The length of 
the tangential spring ^ may increase during the time that the 
condition /* < ^f " is satisfied. The sliding condition is en- 
forced keeping constant the elastic displacement ^ when the 
Coulomb limit condition fjr = /i/^ is reached and fi is the 
interparticle friction coefficient. 

The viscous force is calculated as = —mi'v'^, where 
m is the effective mass of the two particles in contact and 
1/ the damping coefficient. This force takes into account the 
dissipation at the contact and it is necessary to maintain the 
numerical stability of the method. 

A suitable closed set of material parameters for this model 
is the ratio fct/fc„, together with the value of the normal stiff- 
ness kn, the interparticle friction fi, and the ratio ef/e„ be- 
tween the tangential and normal restitution coefficients. 

The random generation of the particles is done by means 
of a Voronoi tessellation using a reference square lattice, 
yielding a set of nearly isotropic polygons. By distorting the 
square lattice in the horizontal and vertical directions, we 
end up with elongated (i.e. anisotropic) particles. The ratio 
between the stretching and contracting factors gives us the 
average aspect ratio A of the polygons, that is used to char- 
acterize the anisotropic shape of the particles. 

In Fig. |2] the different initial sample configurations are 
shown. The isotropic configuration is depicted in Fig. |2^, 
and the anisotropic ones in Fig. for particles stretched in 
the same direction of shearing (horizontal direction, sample 
H) and in Fig. |2j; for particles stretched perpendicularly to 
the shear direction (vertical direction, sample V). 

We use samples of two different sizes, with 256 (16 x 16) 
and 1024 (32 x 32) particles. Periodic boundary conditions 
are imposed in horizontal direction. A constant horizontal 




Fig. 2 (Color online) Samples of (a) isotropic polygons (A = 1) and 
(b) elongated polygons stretched either in the horizontal direction (A = 
2.3 H) or (c) in the vertical direction (A = 2.3 V). 
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Fig. 3 (Color online) Sketch of the shear cell. The system is not al- 
lowed to dilate, i.e. it has fixed h. The sample is sheared using a con- 
stant shear rate 7. Blue (grey) particles induce shear. 

velocity is given to the particles in the top and bottom lay- 
ers so as to impose a constant shear rate 7. These particles 
are not allowed to move in the vertical direction, thus sup- 
pressing the volumetric strain of the system. Furthermore, 
they are not allowed to rotate or move against each other. In 
Fig.|3]a setup of the shear cell is presented for the anisotropic 
sample A = 2.3 H. The shear strain 7 is defined as 

7 - Dx/h, (1) 
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Fig. 4 (Color online) The average kinetic energy in logarithmic scale 
versus the shear strain 7. The stationary value Ko of the kinetic energy 
is obtained from the velocity profile of the paiticles at the steadystate. 
The released energy of the avalanches are calculated using Eq. i|2)- 

cle friction coefficients fi = 0.0, 0.5, 5.0. For simplicity we 
use p = Ig/cm^. 



where is the horizontal displacement of the boundary 
particles and h is the height of the sample, which is kept 
constant. 

In our simple model, polygons represent rocks compos- 
ing the gauge between two tectonic plates and the top and 
bottom boundary particles represent the tectonic plates. We 
start from a perfectly packed configuration in order to rep- 
resent the initial state of the material that is supposed to be 
intact prior to the shear process. 

As described above, the value of the strain rate is of the 
order of 10^^ per year (« 10^^^ s^^). For our numerical 
simulation, this velocity is computationally too expensive. 
For instance, for At ^ Is one needs 10^^ iterations to induce 
a shear strain of 1%. Using a system of 16x16 particles, such 
10^^ iterations would require roughly 1000 years of CPU 
time on a standard P-IV PC. To overcome this, we choose 
a suitable shear rate at which the motion of the system is 
intermittent, i.e. in some regions the system is locked and 
deforms steadily accumulating elastic strain and in others 
the stored energy at the contacts is suddenly released. See 
Fig. ID 

An important issue in this scope is to study the distribu- 
tion of the released energy, spanned over several orders. We 
tested shear rates in the range 10^ - lO"'^ s-i and found that 
a suitable value for the above purpose is 7 = 1.25 x 10~^ 
s^^. Further, to perform the MD simulations using the se- 
lected shear rate, we adjust the parameters of the model in 
order to obtain a time step At requiring a reasonable CPU 
time. Thus, we use fc„ = 400 N/m, e„ = 0.9875, kt/k„ ~ 
1/3, t't/t'n = kt/kn and et/cn = 1.0053, yielding a time 
step of At = 0.005 s. We consider three different interparti- 



3 System response: monitoring avalanches 

The motion of the particles in the interior of the sample is 
not continuous, but has a "stick-slip character". During slip 
a sudden rearrangement of the medium arises as a conse- 
quence of the large relative displacements of the particles. 
We monitor this rearrangement of the system through its 
kinetic energy K. As shown in Fig. |4] the system can be 
in two different states. In the "steady state", K is approx- 
imately equal and less than a low value Ko, shown by the 
horizontal line in Fig. |4l This value Kq is associated with 
the accumulation of elastic strain under the imposed shear. 
If the particles were rigid, we would have Kq = 0, but since 
we are using soft elastic ones, we have Kq > 0. The low 
energy state Kq is punctuated by a series of events where 
kinetic energy rises several orders of magnitude above Kq. 
These are the avalanches. An avalanche begins when K rises 
above Kq, and all subsequent values of K greater than Kq 
are considered to be part of the same avalanche. 

The total released energy Er of one avalanche is the sum 
over the total number N of consecutive values of K above 
the stationary state, namely 

N 

i=i 

where jAt is a proper non-dimensional parameter for in- 
tegrating the kinetic energy (see Fig. |4ll. At the 'stationary 
state' the system is deforming steadily and accumulates en- 
ergy at the particle contacts. This state can be characterized 
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Fig. 5 (Color online) Accumulation of elastic strain and overcome of the strength Ft / Fn of the material prior to the occurrence of an avalanche. 
In (a) the kinetic energy of the system and (b) the ratio Fs/Fn showing the developed strength are shown, with red circles indicating the strain 
value at which the snapshot in (c) and (d) are taken, (c) The configuration and accumulated particle rotation just before the avalanche, (d) The 
elastic strain at the contacts before the avalanche, where the diameter of the red dots is proportional to the value of the deviatoric strain. System 
size: 32 x 32 particles. 



by the value Kq obtained from the average velocity profile 
of the particles at this stage. 

In the case of infinitely rigid particles, subtracting the 
'stationary value' Kq from the kinetic energy of the sys- 
tem, one would obtain zero between successive avalanches. 
Our system is however composed of soft elastic particles, 
and consequently a non-zero value is observed, as shown in 
Fig- HI This non-zero value is a numerical artifact ll24ll stem- 
ming from the calculation of the tangential contact forces, 
the soft elastic nature of the polygons, and the periodic bound- 
ary conditions that can trap some of the energy released dur- 
ing the avalanches. We checked that our present results are 
not affected by this numerical noise and therefore we will 
not consider it in anymore detail. 

The force needed to sustain the constant motion of the 
top and bottom layers can be measured in the simulation. 
In the following, Fg is the shear (horizontal) force applied 
to each wall, and F„ is the normal force. Figure |5] shows 
the occurrence of one avalanches and the associated strain 
accumulation for a system with 32 x 32 particles. We can 
see that the abrupt increment of kinetic energy of the sys- 
tem (Fig.|5t), matches with the fall-off of the strength of the 
material Fg/Fn. 



Figures |5}; and|5}l illustrate two different representations 
of the same sample snapshot, immediately before the ava- 
lanche. Figure |5}; shows the sample configuration and the 
rotation that the particles undergo for a shear band located 
at the center of the sample. The colors of the particles are 
given by their accumulated rotation: the lighter the color the 
bigger the accumulated rotation. Figure |5}l shows the elas- 
tic strain at the contacts, which are represented by red dots 
with a diameter proportional to its strain value. Here, one can 
see that there is a strong localization of elastic strain along 
the shear band. This strain localization weakens the system 
and drives it to failure, since it promotes the occurrence of 
the Coulomb limit condition related to the number of sliding 
contacts (see above). In other words, the weakening of the 
system is due to both the strain localization and the increase 
of the ratio of sliding contacts. 



During the avalanche the system suffers a complete re- 
arrangement in which the old sliding contacts are removed 
from the sliding condition and new contacts are generated. 
This rearrangement marks the beginning of a new stage of 
elastic strain accumulation that drives the system to the next 
avalanche. 
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Fig. 6 (Color online) Log-log plot of the number of avalanches versus 
their released energy Er for (a) different configurations of isotropic 
particles and (b) for different A values. Here /i — 0.5 and the system 
has 16 X 16 particles. Logarithmic binning is used. 



Fig. 7 (Color online) The distribution D(Er) of the released energy 
Er when (a) varying interparticle friction coefficient n with fixed A — 
1 and (b) when varying A with fixed /j. — 0. The system has 16 x 16 
particles and logarithmic binning is used. 



4 The Gutenberg-Richter law in anisotropic granular 
media 

The distribution of earthquake magnitude is described by the 
Gutenberg-Richter law |23l- This law states that the number 
n of earthquakes of magnitude M is proportional to n ~ 
IQ-bM Typically, the value of b is equal to 1 at most places, 
but may vary between 0.8 and 1.5 ll25h . As we will see, 
the exponent b will be an invariant property describing the 
occurrence of avalanches associated with sudden rearrange- 
ments of granular media under very slow shear. 

To this end, we study the possible influence of the for- 
mation and evolution of the shear band on the distribution 
D{Er) of the released energy Ej. during an earthquake. Since 
the magnitude of an earthquake is defined as the logarithm 
of the released energy apart proper constants one finds n ^ 
E^^-, with the exponent c varying typically in the range 0.8 < 

c < 1.1 m. 



In Figure |6^ we show the distribution D{Er) for three 
different initial configurations of isotropic samples, corre- 
sponding to different seeds for the Voronoi Tessellation. All 
the distributions collapse and show a power law behavior 
over almost six orders of magnitude with an exponent of 
c = 0.87 for the fitted straight line, which is within the ob- 
served range of values of the Gutenberg-Richter law. 

In Fig. ISh, the distributions for both isotropic and aniso- 
tropic particles are shown. Similarly, for all samples, the data 
sets are well fitted by a power law with an exponent c rang- 
ing from 0.82 to 0.89, indicating a weak influence of the 
particle shape on the distribution of the released energy. The 
power law holds over six orders of magnitude. 

Similar exponents (0.80 < c < 0.95) are obtained for 
other system sizes in both isotropic and anisotropic cases 
and for the case when one considers the distribution for in- 
dividual particles. From such results, one concludes that in- 
dependent of the anisotropy there is a scale invariance of the 
system response according to the Gutenberg-Richter law. 
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Waiting time t (s) 

Fig. 8 (Color online) Distribution n{t) of waiting times for the se- 
quence of aftershocks in the numerical simulation. Isotropic sample 
A = 1 and anisotropic samples A = 2.3 are presented. The system has 
16 X 16 particles. 

We also study the influence of the friction coefficient fi. 
In Fig.|7^ the distributions for the isotropic samples with dif- 
ferent friction coefficients are plotted. The effect of friction 
in both cases is to slightly increase the exponent c, which 
holds for nearly seven orders of magnitude. The distribu- 
tions for isotropic and anisotropic samples with /^t = are 
presented in Fig. |7J) where no influence of particle shape is 
observed. 



5 Waiting times and Omori's law 

Earthquakes usually occur as part of a sequence of events, 
in which the largest event is called the mainshock and the 
events prior and after the mainshock are foreshocks and af- 
tershocks respectively fl3h . The empirical law that describes 
the behavior of the temporal sequence of avalanches is called 
Omori's law, and states that the number n{t) of aftershocks 
decreases with the inverse of the time interval t spanned 
from the last mainshock as 



where d is an emp irical constant and p varies in the range 
0.7 < p < 1.5 lfl3ll with the most typical values around one. 

Before performing the calculation of waiting times of af- 
tershocks in the system evolution, we have also to precisely 
define 'mainshock' . Our definition is based on empirical ob- 
servations |[l3h . A new event is considered mainshock only 
when its released energy is larger than 1/10 of the released 
energy of the last mainshock. When this happens the se- 
quence of the aftershocks from the previous mainshock is 
considered to be finished and a new sequence is calculated. 

In Fig. |8] the distribution of waiting times for isotropic 
and anisotropic systems are shown. Over more than three 



orders of magnitude all the numerical results can be fitted 
using the expression in Eq. OJ, with exponents p = 1.57 for 
A = 1, p = 1.61 for A = 2.3 V, and p = 0.83 for A = 2.3 H. 
While for A = 2.3 H one observes an exponent within the 
typical range found in fault gouges, for A = 1 and A = 2.3 
V one finds a clear deviation from the observed values. This 
indicates not only that anisotropy should be ubiquitous in 
fault gouges, but also that the most stable configurations are 
the most common, as explained below. 

Using this influence of the initial configuration of aniso- 
tropic samples on the stability of the system, we will next ex- 
plain how to detect at the macro-mechanical level the pres- 
ence of anisotropic particles within the gouge. 

The anisotropic sample A = 2.3 H with particles oriented 
parallel to the shear direction exhibits a more stable configu- 
ration, than the other two cases A = 1 and A = 2.3 V. In this 
sample, the induced torque on the particles is minimized and 
the main deformation modes, sliding and rolling of the parti- 
cles, are highly suppressed by the fixed boundary conditions 
allowing no dilation in vertical direction. The hindrance of 
the deformation modes produces a larger temporal stability 
and also a larger mechanical stability. The larger temporal 
stability makes the occurrence of the events less frequent in 
time, i.e. a slower decay of the waiting times. The larger 
mechanical stability results in a smaller probability of fail- 
ure for a given value of stiffness as discussed in Sec. |6] be- 
low. On the contrary, the configuration of anisotropic sam- 
ples A = 2.3 V, with particles oriented perpendicular to the 
shear direction, maximizes the induced torque on the par- 
ticles and results in a less stable configuration. This con- 
figuration yields smaller temporal and mechanical stability. 
The smaller temporal stability is observed in the decay of 
the waiting times, that is slightly faster than the one of the 
isotropic sample. The smaller mechanical stability of sample 
A = 2.3 V is manifested in the larger probability of failure 
for a given value of stiffness compared to the other samples. 

Therefore, by looking at the decay rate of the aftershock 
sequences one might be able to explain the variation of the 
exponent p in realistic earthquake sequences, by the exis- 
tence of anisotropic gouge in the fault zone. It is important to 
say that for a more realistic representation of the earthquake 
process the crushing of particles should also be taken into 
account. Nevertheless, the absence of particle crushing is a 
suitable approximation for the case of young fault gouges. 



6 Weakening and mechanical stability of the system 

In this section we study the relationship between the occur- 
rence of avalanches and the weakening of the system. The 
weakening process results from the release of energy due to 
previous accumulation of strain at the contact level and con- 
tacts reaching the sliding condition. 

In Fig.|9]we show the relative number of sliding contacts 
Ug, the stored energy at the contacts Eg and the total kinetic 
energy Er. The stored elastic energy Eg at the contacts is 
calculated as Es = 1/2 (fc„(5^ + kt£,'^), where 5 and ^ are 
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Fig. 9 (Color online) Evolution of the relative number Us of sliding contacts, the energy Eg stored at the contacts and the total kinetic energy 
Er as a function of the shear strain 7, for an isotropic sample (A = 1). 



the spring elongations in the normal and tangential direc- 
tions respectively. Figure |9] shows that between avalanches 
the relative number of sliding contacts increases with the 
shear strain. It makes the system weaker and indicates that 
the system is constantly accumulating elastic energy Eg at 
the contacts. The weakening of the system persists until fail- 
ure, where the kinetic energy increases by several orders of 
magnitude. At this stage, the structure of the system is rear- 
ranged, the stored energy at the contacts is released and the 
contacts do not fulfill the sliding condition. All the events in 
the kinetic energy are associated with both drops in the ratio 
Us and drops in the stored energy. 

At the macromechanical level the weakening of the sys- 
tem is observed by looking at the evolution of the shear 
stress with the shear strain. After each stress drop the system 
experiences a rearrangement. This new configuration pro- 
duces a temporal stability, in which the strength builds up. 
At this stage the system sticks and the accumulation of strain 
takes place. In this softening regime the system approaches 
failure and when the strength of the material is overcome, 
the system slips. 

The mechanical stability of the system is studied through 
the stiffness value At/A'^ that the system presents at fail- 
ure. To this end, we calculate the conditional probability 
P{Ae\ At I A^) for the occurrence of an avalanche event 
given a stiffness value At j A^. Since data only provides the 
complementary conditional probability P{At j A^\Ae) we 
use the Bayes theorem 1260, to obtain 



Having a certain time increment dt in the data, the con- 
ditional probability PI^At j A^\Ae) gives the probability of 
observing a stiffness value At j A^ at t — dt when at time 
t an avalanche (symbolized hy Ae) occurs. The conditional 
probability takes into account all the avalanches A e within 
a total time interval of 5 • 10^ s. To this end, we discretize 
the total time interval in time increments of dt = ^7/7- 
Since the shear rate 7 = 1.25 • 10"^ s and A^ = 0.0016, 
the time increment is = 128 s is also constant in our 
simulation. We select the same dt for both isotropic and 
anisotropic systems. The probability P{Ae) is fraction be- 
tween the total number of avalanches and the total number of 
time increments, and similarly P{At/A')) is the probabil- 
ity distribution of At/A^. In Fig. [10] the conditional prob- 
abilities P{AE\AT/A'y) and P{At/A-i\Ae) for both the 
isotropic system A = 1 and the anisotropic system A = 2.3 
are shown. For all the samples, the conditional probability 
P{Ae\At / A'y) decreases logarithmically with the stiffness: 
the stiffer the system the smaller is the probability of failure, 
as shown in Fig II Ob. yielding 



P{AE\AT/A-f) ^ q\og{AT/A-f). 



(5) 



P(AE\AT/A-f) 



P{At/A-i\Ae)P{Ae) 
P{ATlA-i) 



(4) 



Further, from Figs. llOb and 1 10b. anisotropic samples com- 
pared to the isotropic ones explore a different range of stiff- 
ness at failure due to the larger rotational frustration that the 
elongated particles undergo. This is specially true for parti- 
cles oriented along the shear direction (A = 2.3 H), due to 
its more stable structure. The exponents of the tail of the dis- 
tributions are q = -42.1 for A = 1, q = -28.6 for A = 2.3 
V and g = -30.78 for A = 2.3 H (see Fig.HOji). 
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Fig. 10 (Color online) Conditional probability distributions (a) 
P{Ae\At /Aj) of the occurrence of an avalanche Ae given a stiff- 
nesszir/zi7 and (b) P ( At /Aj\AE)oi having a stiffness At /Aj at 
the occurrence of an avalanche. System size 16 x 16. Isotropic A = 1 
and anisotropic samples A — 2.3 are presented. Labels H and V cor- 
respond to horizontal and vertical samples. 



The same analysis can be performed for the occurrence 
of an avalanche Ae given a force ratio Fg/Fn value. In 
Figure [TT] the conditional probabilities P{AE\Fs/Fn) and 
P{Fs/Fn\AE) for the same cases as in Fig. [TOl As shown 
in Fig. [TT^. the conditional probability P{AE\Fs/Fn) in- 
creases approximately linearly with Fs/Fn- The higher the 
mobilized strength Fg / Fn the higher the probability of fail- 
ure. In general, the anisotropic samples are able to mobi- 
lize higher frictional strength than the isotropic sample and, 
for the same force ratio, present a lower probability P{Ae\ 
Fs/Fn). This is certainly due to the influence of particle 
shape anisotropy on the global strength of the material lfl9ll . 

The probability distribution P{Fs/Fn\AE) shown in - 
Fig. [TTb of having a force ratio Fs/Fn at the occurrence 
of an avalanche A e follows typically a normal distribution. 



Fig. 11 (Color online) Conditional probability distributions (a) 
P{Ae\Fs/F„) of the occurrence of an avalanche Ae given a fric- 
tional strength Fs/Fn and (b) P[Fs / Fn\AE) of having a frictional 
strength Fg / F„ at the occurrence of an avalanche. System size 16x16. 
Isotropic A = 1 and anisotropic samples A = 2.3 are presented. 
Labels H and V correspond to horizontal and vertical samples. The 
P[Fs/ Fn\AE) follows a normal distribution, except for sample A = 
2.3 V. Solid lines represent the normal distribution of the data. 



with a mean value larger in the anisotropic samples. Con- 
sidering only the anisotropic cases, while one observes the 
same mean, the standard deviation is larger for the less stable 
configuration V, as expected. 

Finally, in Figure [12] the sliding contact ratio Us and the 
stiffness at failure of the system are presented. The stiffness 
is also compared to the released energy of the avalanches. 
Although no clear correlation between the parameters can 
be observed, the maximum value of stiffness that the sys- 
tem can present is bounded by n^, as illustrated in Fig. 1 12b. 
Larger implies smaller stiffness. In Fig. [T2b one sees a 
weak correlation between the stiffness and the released en- 
ergy, with a slight tendency of Ej. increasing with the stiff- 
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Fig. 12 (Color online) Relationship between (a) sliding contact ratio 
Us VS. stiffness at failure, (b) released energy Er of the avalanche vs 
stiffness at failure. Data correspond to an isotropic sample with size 
16 X 16 particles. 



ness. From these observations, it is clear that the accumula- 
tion of strain at the contacts is not the only important agent 
in the weakening process. Therefore, additional ingredients 
have to be taken into account for a more exhaustive analysis. 



7 Concluding remarks 

In this paper we highlighted the importance of particle shape 
on the occurrence of avalanches in granular systems. To this 
end, we used shear cells with periodic boundary conditions 
to mimic the behavior of tectonic faults with fixed bound- 
aries. 

We found that the dynamics of the granular system is 
characterized by discrete avalanches spannirig several order 
of magnitude similar to crackling noise B \2]\\ . We calcu- 
lated the probability distribution of the energy released in 
avalanches, and found it to be in very good agreement with 



the Gutenberg-Richter law for samples with different parti- 
cle anisotropy and different system sizes. Consequently the 
exponent of the released energy distribution can be seen as 
an invariant property of such systems. 

We also studied the temporal distribution of event se- 
quences after a mainshock. We found that the number of af- 
tershocks decreases with a power of the inverse of time. We 
could fit the sequences of waiting times of the aftershocks 
with the empirical expression and obtained exponents in the 
range 0.7 < p < 1.6, similarly to what is observed in real ob- 
servations according to Omoris law. The anisotropic sample 
H exhibits a larger temporal stability making the temporal 
occurrence of the avalanches sparser, due to its larger frus- 
tration of rotation in the corresponding initial configuration. 

The larger temporal stability observed at the macro-me- 
chanical level can be therefore taken as an indication of the 
existence of anisotropic material within the shear zone. This 
could potentially explain the variation of the exponent p ob- 
served in realistic earthquake sequences. 

The dynamics of the system was also related to the stick- 
slip process ll27l : l28h . When one avalanche begins the system 
slips, while between two successive avalanches, the system 
sticks, accumulates elastic energy, and becomes weaker be- 
cause of the increase of the sliding contacts rig. We charac- 
terized the weakening of the system by the stiffness At/A^ 
and derived the conditional probability P{Ae\At/Aj) for 
an avalanche to occur given a stiffness value At/Aj. We 
found that P{AE\AT/A'y) decreases logarithmically with 
the stiffness and with a decay rate larger for isotropic sam- 
ples than for anisotropic ones. Concerning frictional strength 
we found that the probability of an avalanche to occur in- 
creases with the force ratio Fg/Fn- 

The results concerning the conditional probabilities un- 
covered that anisotropic samples can explore a wider range 
of stiffnesses and force ratios than the isotropic sample. This 
is due to the larger kinematic constraint that anisotropic par- 
ticles undergo during shear. Further, in general, since the ini- 
tial configuration corresponding to sample H is the most sta- 
ble configuration with respect to shearing, it is the one with 
lower probability of failure. 

In previous works concerning the avalanches in granu- 
lar piles 1I29I; [30I1 . some avalanche precursors associated to 
the onset of fluidized regions of sliding contacts were found. 
These fluidized regions were located in the weak network 
of contacts. This weak network comprises contacts trans- 
mitting forces smaller than the average force and therefore 
has a minimal contribution to the stress state Isil; [32ll . It is 
known that [29; 32; 33], while the strong contact network is 
responsible for the strength and stability of the packing, the 
weak contact network plays an important role in the destabi- 
lization process. In this scope, since precursors are expected 
to be related to sharp changes in the micro-structure of the 
granular packing [l34i1 . it would be important to study also 
the network of contact forces as a function of the anisotropy, 
in particular as a function of the major axis orientation. 

Finally, although the results from our numerical model 
show good agreement with the processes observed in na- 
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ture, we are aware of the challenges to have a more realis- 
tic simulation of fault zones [!36]. In future work one should 
also extend this approach to the three-dimensional case, con- 
sider the development of transparent (or absorbing) bound- 
ary condition, since the acoustic emission after an avalanche 
are trapped due to the periodic boundary conditions, and 
also implement the grain fragmentation, since natural earth- 
quakes result from the combined effect of frictional instabil- 
ities and rock fragmentation. 
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